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ABSTRACT 

Temporal changes of X-ray to very-high-energy gamma-ray emissions from 
the pulsar-Be star binary PSR B1259— 63/LS 2883 are studied based on 3-D 
SPH simulations of pulsar wind interaction with Be-disk and wind. We focus 
on the periastron passage of the binary and calculate the variation of the syn- 
chrotron and inverse- Compton emissions using the simulated shock geometry and 
pressure distribution of the pulsar wind. The characteristic double-peaked X-ray 
light curve from observations is reproduced by our simulation under a dense Be 
disk condition (base density ~ 10~ 9 g cm" 3 ). We interpret the pre- and post- 
periastron peaks as being due to a significant increase in the conversion efficiency 
from pulsar spin down power to the shock-accelerated particle energy at orbital 
phases when the pulsar crosses the disk before periastron passage, and when the 
pulsar wind creates a cavity in the disk gas after periastron passage, respectively. 
On the contrary, in the model TeV light curve, which also shows a double peak 
feature, the first peak appears around the periastron phase. The possible effects 
of cooling processes on the TeV light curve are briefly discussed. 

Subject headings: gamma rays:theory-radiation mechanisms :non thermal-stars:Be 
-starts:wind, outfiows-stars:individual (PSR B1259-63) 
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Introduction 



Recent improvement of the techniques of ground-based Cherenkov telescopes has 
increased the number and variety of TeV gamma-ray objects. Five gamma-ray binaries 
have been detected so far, with all of them known to be high mass X-ray binary systems. 
Their common emission mechanism has been vastly investigated since their discoveries. 
Each gamma-ray binary is assumed to be composed of a compact object orbiting around 
a massive star. Among them, PSR B1259-63/LS 2883 is the only system for which the 
compact object has been confirmed to be a pulsar. 



PSR B1259-63 is a 48-ms radio pulsar with a spin down 



power of L sH = 8.2 x 1 35 erg s 1 



1994). However, 



20111 ). 



The spectral type of LS 2883 had been known to be B2V (IJohnston et al. 

a recent precise measurement reports the type to be 09.5V (INegueruela et al. 

This correction implies a change of the assumed distance to the system from 1.5 kpc to 

2.3 ± 0.4 kpc, as well as a change of the effective temperature of the star. The orbit has a 

large eccentricity e of 0.87 and a long period P or b of about 3.4 yr. 



Non-pulsed anc 



2005 



Moldon et al. 



non-thermal emissions from t 



energy ranges (jAharonian et al. 



20111). X-rav (IHiravama et al. 



2005 



re binary in the radio 



1996 ; 



Uchiyama et al. 



(IJohnston et al 



20091 ) and TeV 



20091 ) have been reported, and flare-like GeV 



emissions were detected around th e 2010-2011 periastron passage by the Large Area 



Telescope (LAT) on board Fermi (jTam et al. 



2011 



Abdo et al. 



201lh . The radio pulse 



eclipse of about 5 weeks around the periastron suggests that the pulsar goes to the opposite 
side of the Be-disk plane with respect to the observer, while crossing it plane twice during 



the course. 



he characteristic double-peaked features obs erved in the radio and X-ray 



light curves ( jConnors et al. 



2002; 



Chernyakova et al. 



20061 ) can be mainly attributed to the 



interactions of the pulsar wind and the Be-disk during disk crossings by the pulsar. The 
peak phases or the peak intensities measured, extensively in the radio band in particular, 



-4- 



vary f rom orbit to orbit, tho ugh. The observations for 2004 (jAharonian et al.l 120051 ) and 



2007 (JAharonian et al. 



20091 ) periastron passage indicate that the TeV light curve also 



varies from orbit to orbit. 



The non-thermal emission mechanisms of the system have been studied in 



the framework of leptonic (e.g., 



Sierpowska-Bartosik fc Bednarek 



hadronic models (IKawachi et al. 



Tavani fc Arons 



2008 



1997 



Takata fc Taam 



Khangulvan et al. 



2009 



Dubus et al 



Neronov &: Chernyakova 



2007; 



20101) and 



20071 ). Unfortu- 



20041 ; see also 

nately, Be-disk models which have been adopted so far in this field of research are mostly 
outdated, and are significantly different from the model being widely acc epted by the Be 



star rese arch commun ity c urrently, i.e., the visco u s decr etion disk model 



Lee et al. 



see also 



Porter! (119991 ) and ICarciofi fc Bjorkmanl (120061 )]. Since the shock distance and 



(119911 ): 



its geometry depend on the Be-disk model, it is important to examine the high-energy 
emissions with a more realistic Be-disk model. For instance, disk models with supersonic 
outflows, as adopted in most previous studies, have irrelevantly large ram pressure, and 
thus give rise to false location of the shock. 

Detailed two-dimensional hydrodynamic simulations hav e also been performed to 



Bogovalov et al. 



2008 



2Qll|), 



study the wind- wind collision interaction in this system (e.g., 

providing important information on the shock structure. Since the simulations have been 

limited to 2-D, however, the orbital motion has not been taken into account, which can 

have significant effects near the periastron passage. Moreover, the presence of the Be-disk 

can play an important role to the origin of the high energy emission from this system. In 

such a case, the behavior of the system will inevitably manifest in 3-D and be orbital-phase 

dependent. 



In 



Okazaki et al 



(120 111 hereafter paper I), we have studied the interaction between the 



pulsar and the Be star by, for the first time, carrying out 3-D hydrodynamic simulations 
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using the viscous decretion disk model. We have found that for a Be-disk with typical 
density, the pulsar wind strips off an outer part of the Be-disk, truncating the disk at a 
radius significantly smaller than the pulsar orbit. This prohibits the pulsar from passing 
through the disk around periastron passage, which has been assumed in previous studies. 
In other words, a Be-disk with typical density is dynamically unimportant in this system. 



A large Ha equivalent width of —54 A ( JNegueruela et al.ll201ll ). however, suggests 
that the density of the Be-disk of LS 2883 is much higher than typical. It is, therefore, 
interesting to see how the interaction changes if the Be-disk is much denser and dynamically 
more important than those studied in paper I. Given the double-peaked light curves of PSR 
B 1259-63, a particularly important issue is whether a reasonably high Be-disk density can 
lead to a strong pulsar wind confinement, and hence an enhanced emission, during both 
disk-plane crossings. 

In this paper, which is the second of the series, we study high energy emissions from 
PSR B1259-63/LS 2883 system, based on the results of numerical simulations for different 
values of the Be-disk density. We review our 3-D hydrodynamic simulations in section [5] 
and describe our emission model in section |3j Based on the model light curves, we propose, 
in section 14. 2[ a new interpretation of the observed double-peaked feature of the radiation 
(in particular in the X-ray band). The multi- wavelength emission properties are discussed 
in section 14.31 Finally, we summarize our results in section |5j 

2. Numerical Model for the Hydrodynamic Interaction between the Pulsar 

and the Be Star 



The simulations presented be 
basically identical to that used by 



ow are performed wi th a 3-D SPH code. The code is 



Qkazaki et al. 



( 120021 ) (see also 



Bate et al 



19951 ). except 
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that the current version is adapted to systems with winds and a decretion disk, such as 



PSR B1259-63/LS 2883, and takes into account radiative cooling with the cooling 



gener ated by CLOUDY 90.01 for an optically thin plasma with solar abundances (IFerland 



unction 



19961 ). Using a variable smoothing length, the SPH equations with a standard cubic-spline 
kernel are integrated with an individual time step for each particle. In our code, the 
Be-disk, the Be wind, and the pulsar wind are modeled by ensembles of gas particles of 
different particle masses with negligible self-gravity, while the Be star and the pulsar are 
represented by sink particles with the appropriate gravitational masses. Gas particles which 
fall within a specified accretion radius are accreted by the sink particle. 

To reproduce the Be decretion disk, we inject gas particles just outside of the 
stellar equatorial surface at a constant rate. In paper I, we adopted the injection rate of 
3.5 x lO _9 M yr _1 , which gave rise to a typical disk base density of 10~ 11 g cm' 3 [detailed 
model fit of the observed Ha profile s typically provid es a density between 10 -12 g cm -3 



and several times 10 g cm (e.g. 



Silaj et al. 



20101 )]. In the current study, however, we 



compare the resulting radiation spectrum/light curve for three different injection rates, 
fixing all the other parameters at values adopted in paper I. Particularly, for the mass and 
radius of the Be star, we take values typical for a B2V star to ensure consistency with the 
models studied in paper I. The polar axis of the Be star is tilted from the binary orbital 
axis by 45°, and the azimuth of tilt, i.e. the azimuthal angle of the Be star's polar axis from 
the direction of apastron, is 19°. With this geometry, the pulsar crosses the equatorial plane 
of the Be star at r ~ — 10 d and r ~ +25 d, where r is the orbital phase in days relative 
to the periastron passage. We choose rates of 3.5 x 10 _9 MQyr _1 , 3.5 x lO~ 8 M yr _1 , and 
3.5 x lO _7 M yr~ 1 , corresponding to disk base densities of 10~ 11 g cm -3 , 10~ 10 g cm -3 , 
and 10~ 9 g cm" 3 , respectively. Remarkably, we find that the highest injection rate is 
favored by best reproducing the observed light curve. We note that taking the base density 
of 10 _9 g cm -3 for this system is not unreasonable, given that LS 2883 showed the Ha 
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equivalent width of —54 A in quiescence ( Negueruela et al.ll201ll ). which was one of the 



largest equivalent widths Be stars have ever shown. 

The Be wind and pulsar wind are turned on in the simulation at a certain time after 
the Be-disk has fully developed in the tidal simulation (for more details, see paper I). We 
started the wind simulation at t = 11.44 P or b, 74 days prior to periastron passage. For 
simplicity, we assume that the winds coast without any net external force, assuming in 
effect that gravitational forces are either negligible (i.e. for the pulsar wind) or are cancelled 
by radiative driving terms (i.e. for the Be wind). The relativistic pulsar wind is emulated 
by a non-relativistic flow with a velocity of 10 4 kms _1 and an adjusted mass-loss rate so as 
to provide the same momentum flux as a relativistic flow with the same assumed energy. 
We assume that all the spin down energy L st } = 8.2 x 10 35 erg s _1 goes to the kinetic energy 
of a spherically symmetric pulsar wind. We also assume the Be wind to be spherically 
symmetric with a mass loss rate of 10 -8 M Q yr _1 . 

It is noted that the unshocked pulsar wind is emulated by a non-relativistic flow with 
a momentum flux (L sc ^/c), where c is the speed of light. The global structure of pulsar 
wind under the influence of the Be-wind/disk should be reasonably represented by this 
treatment. This is because it depends primarily on the momentum flux ratio and not 
on whether the pulsar wind is modeled as relativistic (see paper I for details). This will 
justify the application since we estimate the light curve and spectrum by integrating all the 
contributions from each emission region. The resulting light curve and spectrum should not 
depend on the detailed local structure of the shocked pulsar winds. Rather, they should 
depend mainly on the global structure. Thus we conclude the resulting light curve and 
spectrum in this study are robust in spite of our non-relativistic treatment. Of course, 
as stated in Bogovalov et al. (2008, 2011), the relativistic effects can affect the detailed 
structure of unshocked/shocked pulsar wind regions. Thus it is our future plan to extend 



our code to the relativistic regine and examine its influence. The possibility of relativistic 
Doppler boost effect is discussed in section 4.1. 



3. Emission Model for Shocked Pulsar Wind 

We calculate the emission at each orbital phase using data from the 3-D hydrodynamic 
simulations. First, the simulation volume is divided into uniform grids, and at each 
grid point the pulsar wind pressure is calculated by counting only the contribution from 
the pulsar wind particles. Then, from the 3-D distribution of the pulsar wind pressure, 
synchrotron and inverse-Compton emissions are locally evaluated by using an assumption on 
the local magnetic fields and the calculation scheme described in the following subsections. 
In these calculations, we implicitly assume that particles are accelerated through the 
lst-order Fermi mechanism at the pulsar wind termination shock, and have a power- law 
energy distribution. We anticipate that the motion of each SPH particle, which represents 
an ensemble of electrons and positrons, corresponds to the bulk motion of the pulsar wind 
particles represented by the SPH particle. In the down stream region, we assume that the 
motion of the electrons and positrons represented by each SPH particle is randomized and 
the energy distribution is described by the power law function (section 13.2. ip . 

The total emissions are obtained by integrating the local emissions over the whole 
simulation volume. We adopt 100 3 uniform grids and have confirmed that the result does 
not change even if grids with higher resolution are used. Although the integration is taken 
over the whole simulation volume, virtually all the contribution to the high energy emission 
arises from the shocked pulsar wind region, because the unshocked (upstream) region of 
the pulsar wind is too cold to make any significant contribution to the emission. It has 
been suggested that the relativistic bulk motion of the upstream flow may produce a 
considerable high energy emissions via the inverse-Compton process. In this study, however, 
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because our simulations are done in the non-relativistic limit, we do not investigate this 
point. It has been suggested th at the radio emission can come from a larger volume than 



emissions in other wavelengths ( iMoldon et al. 



201ll ). This volume is roughly ~ 100 AU, 



which is about one order of magnitude larger than the size of the present simulation volume 
(~ la ~ 1.5 x 10 14 cm). Estimation of emission in the radio band is hence omitted here and 
postponed to future studies. 

3.1. Magnetic Fields 

Our simulations are performed in the hydrodynamic limit and there are theoretical 
uncertainties involved in determining the magnetic field of the pulsar wind. The magnetic 
field in the shock down-stream is obtained using the Rankine-Hugoniot relations at the 
shock surface and the adiabatic expansion-law of ideal MHD flow as in the cases of steady 
nebulae (e.g. Kennel and Coroniti 1984; Tavani and Arons 1997; Nagataki 2004). The flow 
in our simulations, however, is too dynamic for the above calculation scheme to be applied 
to. Therefore, we adopt another approach frequently used in modeling gamma-ray bursts 
(e.g. Sari et al. 1998; Xu et al. 2011) as follows. On the pulsar side, the total pulsar wind 
pressure P tot at each grid is associated with the local magnetic field B as 

B = y/ri x 87rP tot , (1) 

where t] is fixed to 0.1, which gives a X-ray flux consistent with observations of 
PSR B1259-63/LS 2883 system. 
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3.2. Scheme of Calculation 

3.2.1. Energy distribution of accelerated particles 

The synchrotron and inverse-Compton processes are calculated in the simulation grids 
where the pulsar wind pressure P tot is non-zero and the contribution of the wind particles 
is P g = (1 — 77) -Ptot- We assume that the number of shocked particles per Lorentz factor per 
volume at each grid point is given by a single power law function of 

/(r) = ^r-* (r min < r < r max ), (2) 

where T is the Lorentz factor of the particles. Note that the accelerated particles are 
strictly relativistic in our model, such that their momenta are simply m e cT, where m e is 
the electron rest mass. We will therefore write down their distribution in V instead of 
momentum for convenience. 

The minimum Lorentz factor r m i n is similar to the Lorentz factor of the bulk motion 
of the unshocked flow. At the limit of small magnetization parameter (a <C 1), the latter is 
given by ~ o~lTi, where <tl and Tl are the magnetization parameter and the bulk Lorentz 
factor at the light cylinder radius of the pulsar. With ctl ~ 10 3 and Tl ~ 10 2 " 3 , r min of 
5xl0 5 is obtained. The maximum Lorentz factor r max is determined as r max = min(T s , T c ); 
T s is the Lorentz factor at which the acceleration timescale t a ~ Tm e c/eB with e being 
electron charge is equal to the synchrotron loss timescale t s ~ 9m 3 c 5 /4e B 2 T, whereas 
T c is defined as the Lorentz factor at which the acceleration timescale t a is equal to the 
dynamical timescale of the shocked pulsar wind t c . Here, t c is deduced as ~ /i p /(c/v3) with 
the scale height h p of the gas pressure estimated from the simulation. The typical r max in 
this work is of order 10 8 . 

The energy density e p is related to the gas pressure P„ as e p = 3P g . Using the condition 
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e p = m, 



,c 2 J J Tf(T)dTdQ, we calculate the normalization factor K as 
t ^ 3P f r 2-p~L- P for p ^ 

K = -^^ <{ rmaX " r m in (3) 

[^(rmax/rmin)]" 1 for p = 2 



m P c 



It is well known that the standard lst-order Fermi mechanism in the test-particle 
limit gives p = (r + 2)/(r — 1), where r is the compression ratio at the shock wave. 
Assuming a strong, non-relativistic shock, the compression ratio r = 4 is derived from 
the Rankine-Hugoniot relation, so that the indexp = 2. From previous studies, in general 
we can expect that p > 2 for weak, non-relativistic shocks, and 1.0 < p < 2.2 for strong, 
relativistic shocks, respectively (see, e.g. Longair 1994 and references therein). Limited by 
the spatial resolution, nevertheless, it is very difficult to determine the shock conditions 
in the emission region directly from the hydro simulation. Therefore, unless mentioned 
otherwise, we will adopt p = 2 as our canonical value for the results in this paper. In 
addition, cases for p = 2 ± 0.5 will also be investigated to illustrate the possible effects of 
the uncertainty in p on our results, in particular the multi-wavelength emission spectra. 

3.2.2. Formulae 

The synchrotron power per unit energy emitted by each electron is calculated according 
to Rybicki & Lightman (1979) as 

where E syn = 3her 2 Bsm8 p /4:7rm e c is the typical photon energy and F sy (x) = 

x J Kr,/^(y)dy, where K5/3 is the modified Bessel function of order 5/3. For the pitch angle 

9 P , we use the averaged value corresponding to sin 2 9p = 2/3. The po wer per unit energy 



and p er unit solid angle of the inverse- Compton process is described by (JBegelman &: Sikora 



19870 
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rl P r tic rlrr 1 

<g = v>l (1-^wv^.n., (5) 

where da' /dQ! is the differential Klein-Nishina cross section, f3 = y/T 2 — 1/r, T> = 
r _1 (l — /3cos#i) -1 with 9 1 and 9q being the angles between the direction of the particle 
motion and the propagating direction of the scattered photon and the background photons, 
respectively, h is the Planck constant, Ij, is the background photon field and 9 C expresses the 
angular size of star as seen from the emission point. For the target stellar photons, a soft 
photon field with an effective temperature of 30,000 K is taken, and the Be-disk emission is 
omitted for simplicity. 

For the Earth viewing angle, we assume that the inclination angle of the orbital plane 
with respect to the sky is % ~ 23° and the true anomaly of the direction of Earth is about 
(f) ~ 130° (Johnston et al. 1996; Negueruela et al. 2011). The model spectrum of the 
emission from the shocked wind measured at Earth is calculated as 



P '77 






(6) 



where Ej expresses the summation of the each grid and SVi is the volume of the each grid. 
The optical depth, r 77 for the pair-creation process between the gamma-rays emitted by 
the wind particles and the stellar photons is expressed by 

r 77 = / d£ dE s a 77 dN s /dE s , (7) 

Jo Je c 

where £ is the propagating distance of the 7-ray from the emitted point, dN s /dE s is the 

distribution of the number density of the stellar soft photon and 



3 



<T 77 (E 7 ,-E' S ) = — cr T (l 



16 



V 



2^ 



(3 - w 4 )lni^ - 2t;(2 - v 2 ) 



where a? is the Thomson cross section, v(E 7 ,E s ) = y/l — E c /E 7 and E c = 2(m e c 2 ) 2 /[(l — 
cos^)^] with 77 being the collision angle. We use the distance D = 2.5 kpc in this work. 
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4. Results and Discussion 

In this section, we first summarize the results of the hydrodynamic interaction between 
the pulsar wind and the circumstellar material of the Be star. We then present the 
multi-wavelength light curves and an time-averaged spectrum obtained using the simulation 
data. We will also briefly discuss the cooling processes and comment on our non-relativistic 
emulation of the relativistic pulsar wind. 

4.1. Hydrodynamic Interaction between the Pulsar and the Be star 

As mentioned in section [2j we have run 3-D SPH simulations of hydrodynamic 
interaction in PSR B1259-63/LS 2883 with the base density of the Be decretion disk in 
the range of 10~ 9 — 10~ 11 g cm" 3 . Figures [TJ and [2] show snapshots of the interaction 
between the pulsar wind and the circumstellar material of the Be star in the binary orbital 
plane at two different epochs, 11 days prior to periastron passage (r = —11 d) and 33 
days after it (r = +33 d), respectively. Each figure compares the shock structure in two 
SPH simulations with different disk base densities, po — 10~ 11 g cm -3 (upper panels) and 
Po = 10~ 9 g cm" 3 (lower panels): left and right panels show the distribution of the volume 
density and the pulsar wind pressure, respectively. To clarify the location where the pulsar 
wind is terminated, the right panel also shows the distribution of the volume density in 
contours. These figures highlight the effect of the Be-disk density on the geometry of the 
interaction surface. For po = 10~ n g cm -3 , the pulsar wind easily strips off an outer part 
of the Be-disk, truncating the disk at a radius significantly smaller than the pulsar orbit. 
As a result, the interaction surface is open and covers only a small solid angle around the 
pulsar (Figures [1] and El upper panels), implying that only a small fraction of the pulsar 
wind energy is available for particle acceleration. 
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In contrast, for po = 10~ 9 g cm -3 , the Be-disk has a large enough inertia not to be 
pushed away so easily by the pulsar wind. Thus, as the pulsar approaches the Be-disk 
before periastron passage, the distance between the interaction surface and the pulsar 
rapidly decreases and finally becomes smaller than the scale-height of the disk. The pulsar 
then penetrates the Be-disk, opening a small cavity around it. The disk gas surrounding 
the pulsar terminates the pulsar wind over a large solid angle, converting a large fraction 
of the bulk pulsar wind energy into the energy of shocked particles. After periastron, the 
pulsar approaches the Be-disk again, now moving away from the Be star. Since the pulsar 
wind pushes the disk in the direction along which the disk density decreases, it can move 
the disk gas more easily than it could before periastron. As a result, a slowly expanding 
shell is formed around the pulsar, which terminates the pulsar wind again over a large solid 
angle, efficiently converting the pulsar wind energy into the energy of shocked particles 
(Figure |2l lower panels). The details of these hydrodynamic simulations will be discussed 
in a subsequent paper (Okazaki et al. 2012, in prep) 

4.2. X-ray Light Curves 

Figure E] presents the calculated light curves of the X-ray flux (1-10 keV energy bands) 
together with the observed fluxes taken from Neronov& Chernyakova (2007) as a function of 
days relative to the periastron passage, r. The solid, dashed and dotted lines are the results 
for disk base densities of po = 10~ n g cm -3 , 10~ 10 g cm" 3 and 10~ 9 g cm -3 , respectively, 
with the typical value of the power index p =2 for the accelerated particles. 

During the pre-periastron period, up to r = — 11 d, the flux does not depend on the 
disk base density. This is because the pulsar wind interacts mainly with the stellar wind, 
for which we assume an identical mass loss rate of lO _8 M yr _1 . The epochs during 
which the pulsar crosses the Be-disk plane are estimated by the simulation to be about 11 
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days prior to the periastron passage (r ~ —lid) and 25 days after it (r ~ +25 d). The 
interaction between the pulsar wind and Be-disk is expected to be the strongest during the 
disk crossings. 

Remarkably, however, the X-ray flux for cases of disks with typical base densities 
Po = 10~ 10 — 10~ n g cm -3 does not peak at the timing of disk crossings, but shows a 
maximum intensity at the periastron. On the other hand, the X-ray flux for the highest 
density case (po = 10~ 9 g cm -3 ) increases distinctively around the phases of the disk 
crossings, resulting in a double-peaked structure. 

The dependence of the light curves on the Be-disk density reflects that of the shock 
geometry and resultant conversion efficiency from the spin down power to the internal 
energy of the shocked pulsar wind (Figured]). For a typical disk density, we see that the 
pulsar wind can easily truncate the Be-disk at a radius smaller than the pulsar's orbit. 
Therefore, the solid angle as measured from the pulsar, over which the pulsar wind is 
stopped by the Be-disk, is small. In this case, the intensity of synchrotron emission is 
highest at the periastron because the magnetic fields of the shocked pulsar wind are highest 
there (see equation (1)). 

The pulsar wind cannot dismiss the disk with the higher density (po = 10~ 9 g cm" 3 ), 
i.e., with a large inertia. The shock is pushed back toward the pulsar in the first disk plane 
crossing (r ~ —lid). After the periastron, the pulsar approaches the Be-disk again, but 
now moving away from the Be star. Since the pulsar wind pushes the disk in the direction 
along which the disk density decreases, it displaces the disk gas more easily than it did 
in the pre-periastron crossing. As a result, a slowly expanding shell is formed around 
the pulsar, which terminates the pulsar wind over a larger solid angle (r ~ +33 d). The 
particles at the shock obtain energy from the bulk kinetic energy of the unshocked pulsar 
wind, which is in turn provided by the pulsar spin down power. At these disk-plane crossing 
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phases, the conversion efficiency from the pulsar spin down power to the shock-accelerated 
particle energy drastically increases, which leads to the increase in the X-ray flux, most 
notably in the post-periastron crossing phase (Figure H]). 

We would like to stress that no fine-tuning has been done in our study for the sake of 
reproducing the observed X-ray light curve. Thus it is remarkable that the calculated flux 
level and double-peak phases for the case of po = 10~ 9 g cm" 3 and r\ = 0.1 in equation (JT]) 
turn out to be very similar to those from observations, as shown in Figure [31 Our model 
sheds new light on the disk density as a probe of the high energy emission mechanisms. 
Interestingly, the favored density in this study is higher than the typical one for Be-disks. 

Our model possesses a few profound features which distinguish it from previous studies. 
First, the present simulations predict that the conversion efficiency varies with the orbital 
phase, while the previous studies have assumed a constant con version efficiency over the 



whole orbit (e.g. 



Tavani &: Arons 



1997 ; 



Takata fc Taam 



20091 ). Second, in the present 
model, the conversion efficiency and therefore the X- ray flux acquire t he m aximum values 



when the pulsar wind interacts with the Be-disk. In 



Tavani fc Aronsl (119971 ). on the other 



hand, the inverse-Compton cooling process which dominates over other co oling; processes is 



crucia l to r eproduce th e obser ved decrease of the flux near the periastron. 



Takata &: Taam 



f l2009h and 



Kong et al. 



( I201ll ) invoked the model that the pulsar wind parameters at the 



shock (e.g. the a parameter) vary throughout the orbital phase. 



4.3. GeV/TeV Light Curves 

Figure |5] shows the calculated light curves for 0.1-100 GeV and TeV (>300GeV) energy 
bands. We find that the gamma-ray light curves show the double-peaked structures as 
in their X-ray counterpart. For 0.1-100 GeV light curves, the peaks align with those in 
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the X-ray bands, because both X-ray and 100 MeV-1 GeV emissions are produced by 



the s ynchrotron process. In the recent results from Fermi fjTam et al. 



2011 



Abdoetal 



201ll ). however, it is likely that the phases of the peaks in the X-ray and the GeV 
bands do not align with each other. Furthermore, Fermi has not detected a strong 
and sharp peak before the periastron, which is seen in the present model. Hence, 
in the present framework, there are some discrepancies between the properties of the 
calculated and observed GeV light curves. We would like to note that the 0.1-1 GeV 
flux particularly depend on the magnetic field, because the synchrotron spectrum from 
the shocked particles has a cut-off around 1-200 MeV. In the present model, if the 
maximum Lorentz factor is determined by the balance between the synchrotron cooling 
timescale and the acceleration timescale, the maximum Lorentz factor of the particles is 
expressed as T s ~ {Qm? e c A / Ae^ B) 1 / 2 , which indicates the synchrotron photon energy of 
E syn = 3T 2 eBsm6 p /(47rm e c) ~ 3 5/2 m e c 3 /i/(2 7/2 vre 2 ) ~ 200 MeV. For a lower magnetic field, 
the maximum Lorentz factor may be determined by the balance between the acceleration 
timescale and the dynamical timescale of the shocked pulsar wind. In such a case, the cut-off 
energy of the synchrotron radiation is E < E syn ~ 200 MeV, and the 0.1-1 GeV flux is 
sensitive to the magnetic field. Furthermore, contributio ns by emissions from the high-order 



generated pairs (ISierpowska-Bartosik fc Bednarek 



from the unshocked pulsar wind (e.g., 



2008T) and th e inverse- Compton emission 



Khangulyan et al. 



2011al ) to the observed spectrum 



below the TeV band have been pointed out. Therefore, a more detailed modeling of the 
0.1-100 GeV emission process is necessary for a closer comparison with the Fermi results. 

For the TeV emissions, which are produced by the inverse-Compton process, we can 
see in Figure [5] that the second peak of the light curve aligns with that of the X-ray band, 
whereas the first peak comes around near the periastron. The timing of the first peak of the 
TeV lightcurve reflects the fact that the stellar soft-photon field strength at the emission 
region (that is, the shock) reaches its maximum during periastron passage. As we described 
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in section \A.2\ though the conversion efficiency from the spin down energy to the particles 
energy of the shocked wind decreases as the pulsar moves toward the periastron, the effect 
of the increasing strength of the soft photon fields toward the periastron outweighs the 
decrease of the conversion efficiency. The periastron peak does not well reproduce the 
observations, at least of the 2004 periastron passage. 

It is shown in Figure 5 that the light curve in the > 300 GeV band suffers from 
photon-photon absorption where soft photons are assumed to come only from the Be star. 
In the pr esent calculation, we ha ve taken a spherically symmetric stellar photon field. 
However, 



Negueruela et al. 



( 1201 ll ) estimated a stellar temperature of T c g ~ 3.4 x 10 4 K 
at the pole and T e g ~ 2.75 x 10 4 K on the equator, which implies that the stellar photon 
field does depend on the latitude. Such a latitudinal dependence of the stellar photon 
field could affect the peak phase of > 300 GeV emissions. As mentioned in section [31 
the current model also omits the contribution from the disk emission, but it can become 
important for the spec t rum a nd the absorption of the GeV-TeV photons. For example, 



van Soelen fc Meintjed (120111 ) studied the effects of the IR excess on the spectrum and 
found that the GeV gamma-ray flux can increase by a factor > 2. When we include the 
contribution fro m the Be-disk emission, the ab sorption effect could be more than doubled, 



as inferred from 



van Soelen fc Meintjesl (1201 ll ). Due to this effect, the flux of > 300 GeV 



gamma-rays around the periastron phase may be suppressed. As a result, the peak phase of 
the > 300 GeV gamma-ray light curve may be shifted from the periastron phase to a phase 
prior to the periastron. We are planning to investigate this effect using a Monte-Carlo 
radiative transfer code as our next step. 
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4.4. Multi-wavelength Spectra 

In Figure [61 we show the spectra during the periastron passage in mult i- wavelength. 
On top of them, we plot the observed spectra measured by various instruments. To account 
for the systematic uncertainty in our emission model on the shock conditions and hence 
the energy distribution of the underlying accelerated particles, two model spectra are 
calculated using p=1.5 (left panel) and p=2.5 (right panel), respectively. In the figures, the 
solid, dashed and dotted lines represent the spectra averaged over three different period 

r = -40 d 10 d, r = -10 d - +20 d, and r = +20 d - +50 d, respectively. Because 

we stopped our simulation at r = +50 d, we cannot calculate the spectrum averaged over 
the period covered by Fermi observation. The model spectra below and above 1 GeV 
correspond to the synchrotron emission and the inverse-Compton process, respectively. In 
Figure [6j we can see a typical flux and spectral shape for each of different power indexes. 
The spectral slope does not change much from phase to phase if the power index is constant, 
although the flux level varies. Thus having model spectra will help us to infer the actual 
power index of the energy distribution of the shocked wind particles from the photon index 
of observed spectra. 

The present model shows a spectral break of the synchrotron radiation around 
1-10 keV, which corresponds to the minimum Lorentz factor of the particles, r m j n = 5 x 10 5 . 



This spectra 



feature is consistent with the break around 1 keV measured by the SUZAKU 



observation ( Uchiyama et al. 



120091 ) . As discussed in section 4.3 and also seen in Figure El the 
present model expects that the synchrotron spectrum extends up to the maximum photon 
energy of E syn = 27m e c 3 h/ (16ire 2 ) ~ 200 MeV, which corresponds to the Lorentz factor 
of the particles, T max ~ (§m 2 e c 4 j^Bf' 2 = 3.6 x 10 8 (7y/0.1)- 1 / 4 (P tot /0.1dyne cm" 2 )- 1 / 4 , 
where we used the equation of (JT]). 

Although the flare-like GeV emissions detected by the Fermi telescope may not be 
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compared with the result of the present simple calculation, where we ignore the effects of 
the radiative cooling. We would like to remark that a smaller power index of the particle 
distribution is preferred to explain the observed flux of >100 MeV emissions. In fact, 
Takata & Taam (2009), who fit the observed X-ray flux and the photon index for various 
orbital phases, pointed out that the expected >100 MeV flux for the smaller power index 
can be higher than the Fermi sensitivity. A more detailed modeling for the spectrum and 
light curve in GeV energy bands will be done in our subsequent studies. 

In the present calculation, the inverse-Compton process between the accelerated 
particles with T min > 5 x 10 5 and the stellar photons of energy ~ 1 eV takes place in the 
Klein-Nishina regime. In such a case, the inverse-Compton spectrum calculated with the 
particle power index p is given by EFe oc E l ~ p above the energy ~ r min m e c 2 ~2x 10 11 eV. 
As we can see in Figure El the observed spectrum EFe oc E~ x by H.E.S.S. (filled circles 
and triangles) indicates the power low index p of the distribution of the scattering particles 
is p ~ 2. The present model also predicts that there is a change of the spectral slope at 
~ 5 x 10 11 eV. 

4.5. Effect of Radiative Cooling 

In the shocked regions with strong magnetic fields, the synchrotron cooling timescale 
t s of the accelearated particles may be less than the crossing time t c of those regions. In 
the present framework, using equations (fllh the total pressure and hence the magnetic 
field strength become largest when the pulsar penetrates the disk, so that the effect of 
synchrotron cooling may not be neglected at that phase. In addition, the inverse-Compton 
cooling time may be less than t c at phases near the periastron. These cooling processes 
can affect the resulting spectrum and light curve. For example, if the inverse-Compton 
process dominates the other cooling processes of the TeV particles, the X-ray emission via 
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synchrotron radiation is weaker in denser soft photon fields. On the other hand, if the 
synchrotron process is the dominant cooling process of the TeV particles, the TeV emissions 
via the inverse-Compton process is weaker in regions with stronger magnetic fields. 

For PSR B1259-63/LS 2883, however, we expect that the cooling processes have no 
important effect on the X-ray light curve. The synchrotron radiation occurs in the slow 
cooling regime if the particle's Lorentz factor is smaller than r sy = (9m 3 c 6 )/(4e 4 I? 2 £), where 
£ is the size of the emission cavity. With B ~ 0.1 Gauss and I ~ 10 13 cm as typical values 
near the periastron, the critical energy below which the synchrotron cooling process can 
be ignored is E ~ 7.62 x 10 7 (i?/0.1 G)~ 3 (l/10 13 cm)~ 2 ) eV, w hich is way above the X-ray 



energy band. Moreover, according to 



Tavani &: Arond (119971 ). the inverse-Compton process 



between the accelerated particles and stellar photons enhances the double-peaked structure 
in the X-ray light curve. Hence, our model X-ray light curves remain robust even if the 
detailed cooling processes will be taken into account. 

5. Summary 



In our previous study (jOkazaki et al 



201ll ). we have developed 3-D SPH simulations of 



the interaction between the pulsar wind and the Be-disk and wind in the gamma-ray binary 
PSR B1259-63/LS 2883. In this paper, we investigated the high-energy emissions from 
the shocked pulsar wind, calculating the synchrotron radiation and the inverse-Compton 
process on the basis of the simulated shock geometry and pressure distribution of the pulsar 
wind. The current study revealed that the observed double-peaked X-ray light curves are 
reproduced only if the Be-disk is denser than typical (with base density ~ 10 _9 g cm -3 ). 
The pre- and post-periastron X-ray peaks appear respectively when the pulsar passes 
through the disk prior to the periastron, and when the pulsar wind creates a cavity in the 
disk gas after the periastron, in both cases terminating the pulsar wind over a large solid 
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angle around the pulsar. 

On the other hand, in the model TeV light curve, which also shows a double 
peak feature, the first peak appears around the periastron, which will disagree with 
the 2004 H.E.S.S. observation showing the first peak located at a pre-periastron phase. 
In a subsequent paper, we will study whether the effects of the disk emission to the 
inverse-Compton process and photon-photon absorption process can shift the first peak to 
a phase prior to the periastron passage. 
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Fig. 1. — Snapshots taken from two SPH simulations for different base densities of the Be- 
disk: po = 10~ n g cm -3 (upper panels) and po = 10~ 9 g cm -3 (lower panels). The left and 
right panels show the volume density and the pulsar wind pressure, respectively, at r = — 11 d 
(11 days prior to the periastron). In the left panels, the locations of the Be star (the central 
white dot) and the pulsar (the white dot above the Be star) are also shown, while in the 
right panel, the volume density is denoted by contours. The density varies by one order of 
magnitude per contour. In each panel, N\, N2, and N$ annotated at the lower left corner 
are the numbers of particles in the Be wind, the pulsar wind, and the Be-disk, respectively. 
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Fig. 2. — Same as Figured! but at r = +33 d (33 days after the periastron) 
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Fig. 3. — The light curves in the 1-10 keV energy band. The solid, dashed and dotted 
lines are the results for base densities of po = 10" 11 gem" 3 , 10" 10 gem" 3 and 10 -9 gem" 3 , 
respectively. We assumed a power index p = 2 for the accelerated particles. Overlaid in the 
plot are data points from Neronov & Chernyakova (2007). 
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Fig. 4. — Energy budget for the radiation vs. the orbital phase. The vertical line corresponds 
to the physical quantity Y 1 i(c/3)P to t,i5V i , where 5Vi is the volume of the grids. The solid, 
dashed and dotted lines are results for base densities of p — 10~ n gcm~ 3 , 10^ 10 gem -3 and 
10 -9 gcm~ 3 , respectively. 
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Fig. 5. — The model light curves for different energy bands for po = 10~ 9 g cm -3 . The upper 
and lower panels show the light curves in the 0.1-100 GeV and >300 GeV energy bands, 
respectively. In the lower panel, the dotted line represents the light curve without accounting 
for the effect of photon-photon absorption. The results are for power index p = 2 for the 
accelerated particles. 
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Fig. 6. — Spectrum of the emissions during the periastron passage. The left and right panels 
show the model spectra for p = 1.5 and p = 2.5, respectively. The solid, dashed and dotted 

lines represents the model spectra averaged over three different periods r = —40 d 10 d, 

t = —10 d — h20 d, and r = +20 d — 1-50 d, respectively. The data points are taken from 
Aharonian et al. (2005, 2009) for the H.E.S.S. observations between February and May 2004 
and between April and August 2007, and from Abdo et al. (2011) for the Fermi observations 
and Swift observation of the 2010-2011 periastron passage. 



